prepTexL <- function(mod){
    
    if(class(mod) == "glm") {print("Use clustered std errors")
      } else if(class(mod) == "polr"){
    
    out_mod <- lmtest::coeftest(mod, vcov = sandwich::vcovCL(mod, cluster = df_L$dst))
    
    coef_names <- names(out_mod[,1])
    est <- out_mod[,1]
    se <- out_mod[,2]
    pval <- out_mod[,4]
    aic <- stats::extractAIC(mod)[2]
    deviance <- mod$deviance
    N <- mod$n
    
    gof <- c(aic, deviance, N)
    gof_names <- c("AIC", "Deviance", "N")
    gof_decimal <- c(T,T,F)
    
    texreg::createTexreg(coef.names = coef_names,
                         coef = est,
                         se = se,
                         pvalues = pval,
                         gof = gof,
                         gof.names = gof_names,
                         gof.decimal = gof_decimal)
    } else if(class(mod) == "glm.cluster"){
    
    coef_names <- names(mod$glm_res[["coefficients"]])
    est <- mod$glm_res[["coefficients"]]
    se <- sqrt(diag(mod$vcov))
    pval <- 2*pnorm(abs(est/se), mean = 0, sd = 1, lower.tail = F)
    aic <- mod$glm_res$aic
    deviance <- mod$glm_res$deviance
    N <- length(mod$glm_res$fitted.values)
    
    gof <- c(aic, deviance, N)
    gof_names <- c("AIC", "Deviance", "N")
    gof_decimal <- c(T,T,F)
    
    texreg::createTexreg(coef.names = coef_names,
                         coef = est,
                         se = se,
                         pvalues = pval,
                         gof = gof,
                         gof.names = gof_names,
                         gof.decimal = gof_decimal)
  } else {print("Unknown model")}

}


prepTexH <- function(mod){
  
  if(class(mod) == "glm") {print("Use clustered std errors")
  } else if(class(mod) == "polr"){
    
    out_mod <- lmtest::coeftest(mod, vcov = sandwich::vcovCL(mod, cluster = df_H$dst))
    
    coef_names <- names(out_mod[,1])
    est <- out_mod[,1]
    se <- out_mod[,2]
    pval <- out_mod[,4]
    aic <- stats::extractAIC(mod)[2]
    deviance <- mod$deviance
    N <- mod$n
    
    gof <- c(aic, deviance, N)
    gof_names <- c("AIC", "Deviance", "N")
    gof_decimal <- c(T,T,F)
    
    texreg::createTexreg(coef.names = coef_names,
                         coef = est,
                         se = se,
                         pvalues = pval,
                         gof = gof,
                         gof.names = gof_names,
                         gof.decimal = gof_decimal)
  } else if(class(mod) == "glm.cluster"){
    
    coef_names <- names(mod$glm_res[["coefficients"]])
    est <- mod$glm_res[["coefficients"]]
    se <- sqrt(diag(mod$vcov))
    pval <- 2*pnorm(abs(est/se), mean = 0, sd = 1, lower.tail = F)
    aic <- mod$glm_res$aic
    deviance <- mod$glm_res$deviance
    N <- length(mod$glm_res$fitted.values)
    
    gof <- c(aic, deviance, N)
    gof_names <- c("AIC", "Deviance", "N")
    gof_decimal <- c(T,T,F)
    
    texreg::createTexreg(coef.names = coef_names,
                         coef = est,
                         se = se,
                         pvalues = pval,
                         gof = gof,
                         gof.names = gof_names,
                         gof.decimal = gof_decimal)
  } else {print("Unknown model")}
  
}
